! created 08/30/2006 by Roger Marchand
  
  subroutine histogram_radar_reflectivities( &
		n_total_columns,	&	! number of profiles
		n_lev,			&  	! number of vertical layers in dBZe field
		dBZe,			&	! Radar Reflectivity

		zbound, 		&	! location of z-grid boundaries in reflectivity field
 
		n_dBZe_boundary_pts, 	&	! define dBZ axis to use in histogram
		dBZe_min,		&	! center value given by
		d_dBZe,			&	!   dBZe_boundaries = dBZe_min+i*d_dBZe,i=0,n_dBZe_boundary_pts-1
		
		n_radar_grid_pts,	&	! define Height axis used in historam
		radar_fixed_grid,	&	! desired height grid for histogram
						! points should be fixed distance apart.

		radar_z_res,		&	! resolution of radar volume centered at
					 	! the point in radar_fixed_grid.
						! NOTE: this value may not be the same as
						! spacing in the radar_fixed grid.

		hist_height_vs_dBZe	&	! output
  )

  implicit none

  real*8, dimension(n_total_columns,n_lev) :: dBZe			

  real, dimension(n_lev+1) :: zbound		

  real, dimension(n_radar_grid_pts) :: radar_fixed_grid  

  integer ::	n_radar_grid_pts, &
		n_dBZe_boundary_pts, &

		n_total_columns, &	 ! size of input Radar Refectivity field
		n_lev 
				
  real ::		radar_z_res, &	! size of resolution volume (not necessarily same as samplin in 
			dBZe_min, &
			d_dBZe	
  !
  ! Output
  !
  real, dimension(n_radar_grid_pts,n_dBZe_boundary_pts-1) :: hist_height_vs_dBZe

  ! local variables
  real :: 	layer_weight(n_lev), &
		layer_dBZe(n_total_columns), &
		sum_Ze(n_total_columns)

  integer	i,m,model_level_ind, &
		dBZe_index, &
		lower_bound_level,upper_bound_level

	! the following codes assume zbound Decrease monontonically
	! i.e. zbound(n_lev+1) is closest bin to the ground	
	! flip if needed
	if(zbound(1)<zbound(2)) then
		zbound=zbound(n_lev+1:-1:1);
		dBZe=dBZe(n_lev+1:-1:1,:);		
	endif

	! the following codes assumes the height points INcrease monontonically
	! i.e. radar_fixed_grid(1) is closest bin to the ground	
	if(radar_fixed_grid(2)<radar_fixed_grid(1)) then
		radar_fixed_grid=radar_fixed_grid(n_radar_grid_pts:-1:1);
	endif
	
	hist_height_vs_dBZe = 0		! initialize histogram

	upper_bound_level=n_lev+1	! lowest point in model grid
	lower_bound_level=n_lev+1	

	do i=1,n_radar_grid_pts

		! find location of top or radar resolution volume boundaries
	
		! decrease model level index until model z boundary is higher than the radar grid
		do while (radar_fixed_grid(i)+0.5*radar_z_res>zbound(upper_bound_level) )
		
			upper_bound_level=upper_bound_level-1

			if(upper_bound_level<1) then
				print *,'Radar grid higher than model grid !'
				stop
			endif
		enddo

		! decrease model level index until model lower z boundary is HIGER than the radar grid
		do while (radar_fixed_grid(i)-0.5*radar_z_res>zbound(lower_bound_level) )

			lower_bound_level=lower_bound_level-1
		enddo

		! generate sub_array of dBZe for the fixed radar level
		if( lower_bound_level == upper_bound_level) then

			! then the radar level is entirely with in one model level

			if(lower_bound_level>n_lev) then
			
				! the entire model resolution volume is below the model lowest layer 
				layer_dBZe=-999
			else
				! put this layer in the histogram ....
				layer_dBZe(:)=dBZe(:,upper_bound_level)
			endif
		else
			! need to weight two or more model levels
			
			! write(*,*) upper_bound_level,radar_fixed_grid(i)+0.5*radar_z_res,zbound(upper_bound_level+1)
			! write(*,*) lower_bound_level,radar_fixed_grid(i)-0.5*radar_z_res,zbound(lower_bound_level)

			layer_weight(upper_bound_level) = &
				( (radar_fixed_grid(i)+0.5*radar_z_res) - &
				  zbound(upper_bound_level+1) ) &
				/ radar_z_res
				
			do model_level_ind=upper_bound_level+1,lower_bound_level-1

				layer_weight(model_level_ind) = &
					( zbound(model_level_ind) - &
				   	  zbound(model_level_ind+1) ) &
					/ radar_z_res 
			enddo

			if(lower_bound_level>n_lev) then

	   		   ! radar volue extends below model ... add weight to lowest model bin
			   layer_weight(lower_bound_level-1) = &
				layer_weight(lower_bound_level-1) + &
				( zbound(lower_bound_level) - &
				   (radar_fixed_grid(i)-0.5*radar_z_res) ) &
				/ radar_z_res
			else
			   layer_weight(lower_bound_level) = &
				( zbound(lower_bound_level) - &
				   (radar_fixed_grid(i)-0.5*radar_z_res) ) &
				/ radar_z_res
			endif

			! write(*,*) 'level =', i, 'we=', layer_weight(upper_bound_level:lower_bound_level)
			
			! sum Radar power (Ze)
			sum_Ze=0
			m=0
			do model_level_ind=upper_bound_level,lower_bound_level
 
			   if(model_level_ind<=n_lev) then
	
				m=m+1

				! convert dBZe to Ze and sum weighted contribution
				sum_Ze=sum_Ze+10**(dBZe(:,model_level_ind)/10) &
					*layer_weight(model_level_ind)
			   endif

			enddo

			layer_dBZe=10*log10(sum_Ze/m+1E-30)
					
		endif

		! write(8,*) layer_dBZe

		! add dBZe for this layer to the histogram
		do m=1,n_total_columns

			if(layer_dBZe(m) >= dBZe_min .and. layer_dBZe(m)<1000) then

				dBZe_index=1+(layer_dBZe(m)-dBZe_min)/d_dBZe

				if(dBZe_index.gt.n_dBZe_boundary_pts-1) then
					dBZe_index=n_dBZe_boundary_pts-1
				endif
	
				hist_height_vs_dBZe(i,dBZe_index) =  &
					hist_height_vs_dBZe(i,dBZe_index) + 1	
			endif
		enddo
	enddo
	
	hist_height_vs_dBZe=hist_height_vs_dBZe/n_total_columns		

  end subroutine histogram_radar_reflectivities
